//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __utils_hpp__
#define __utils_hpp__

#include <vector>
#include <list>
#include <stdio.h>
#include <string.h>

template <class T>
class array16
{
public:
  array16 () : mem_ptr(0), mem_begin(0), mem_end(0), len(0) { }

  void resize (size_t s)
  {
    unsigned char *old_mem_ptr = mem_ptr;
    T *old_mem_begin = mem_begin;
    size_t old_len = len;

    len = s;
    if ( old_len < len )
      {
        mem_ptr = (unsigned char *)malloc (16 + (len + 1) * sizeof (T));
        size_t n = ((size_t)mem_ptr)%16;
        mem_begin = (T *)(mem_ptr + n);
        mem_end = mem_begin + len;
        if ( old_mem_ptr )
          {
            memcpy (mem_begin, old_mem_begin, (len < old_len ? len : old_len) * sizeof (T));
            free (old_mem_ptr);
          }
      }
  }

  void erase (size_t index)
  {
    memcpy (mem_begin + index, mem_begin + index + 1, (len - index - 1) * sizeof (T));
    resize (len - 1);
  }

  void release ()
  {
    if ( mem_ptr )
        free (mem_ptr);
    mem_ptr = 0;
    mem_begin = mem_end = 0;
    len = 0;
  }

  T& operator[] (size_t i) { return mem_begin[i]; }
  const T& operator[] (size_t i) const { return mem_begin[i]; }

  T* begin() { return mem_begin; }
  T* end() { return mem_end; }
  const T* begin() const { return mem_begin; }
  const T* end() const { return mem_end; }

  size_t size () const { return len; }

  ~array16 ()
  {
    release ();
  }

private:
  unsigned char *mem_ptr;
  T *mem_begin;
  T *mem_end;
  size_t len;

  array16 (const array16&); // we are a dumb class, go away
  array16& operator= (const array16&);
};

template <class T, size_t CI = 32>
class allocator
{
public:
  allocator() : it(0), end(0), memsize(0) { }

  ~allocator ()
  {
    clear();
  }

  void clear ()
  {
    while ( !cache.empty() )
      {
        free (cache.front());
        cache.pop_front();
      }
    it = end = 0;
    memsize = 0;
  }

  void reserve (size_t n)
  {
    if ( (size_t)(end - it) < n )
      {
        size_t bc = sizeof (T) * n;
        memsize += bc;
        it = (T *)malloc (bc);
        end = it + n;
        cache.push_back (it);
      }
  }

  T * get()
  {
    if ( it == end )
        reserve (CI);
    return it++;
  }

  T * get(size_t n)
  {
    reserve (n);
    T *ret = it;
    it += n;
    return ret;
  }

  T * getz()
  {
    T *ret = get();
    memset (ret, 0, sizeof(T));
    return ret;
  }

  T * getz(size_t n)
  {
    reserve (n);
    T *ret = it;
    it += n;
    memset (ret, 0, sizeof(T) * n);
    return ret;
  }

  size_t size() const { return memsize; }

private:
  typedef std::list<T *> cache_t;
  cache_t cache;
  T *it;
  T *end;
  size_t memsize;
};

class FPWriter
{
public:
  FPWriter (FILE *fp) : fp(fp) { }

  void write (const void *ptr, size_t len)
  {
    fwrite (ptr, len, 1, fp);
  }

  void write_char (char c)
  {
    fwrite (&c, 1, 1, fp);
  }

  void write (size_t v)
  {
    unsigned int n = v;
    fwrite (&n, sizeof (n), 1, fp);
  }

  void write (float v)
  {
    fwrite (&v, sizeof (v), 1, fp);
  }

  void write (long double v) { write ((double)v); }

  void write (double v)
  {
    fwrite (&v, sizeof (v), 1, fp);
  }

  void write (const vect3d& v)
  {
    write (prec2double (v.x));
    write (prec2double (v.y));
    write (prec2double (v.z));
  }

  ~FPWriter () { if ( fp ) fclose (fp); }

private:
  FILE *fp;
};

class FPReader
{
public:
  FPReader (FILE *fp) : fp(fp), error(false) { }

  void read (void *buf, size_t len)
  {
    do_fread (buf, len, 1);
  }

  void read (double& v)
  {
    do_fread (&v, sizeof(v), 1);
  }

  void read (long double& v)
  {
    double dv;
    read (dv);
    v = dv;
  }

  void read (float& v)
  {
    do_fread (&v, sizeof(v), 1);
  }

  void read (vect3d& v)
  {
    double x, y, z;
    read (x);
    read (y);
    read (z);
    v.x = x;
    v.y = y;
    v.z = z;
  }

  void read (bool& v)
  {
    size_t bv;
    read (bv);
    v = bv != 0;
  }

  void read (size_t& v)
  {
    unsigned int n;
    do_fread (&n, sizeof(n), 1);
    v = n;
  }

  void read_char (char& c)
  {
    do_fread (&c, 1, 1);
  }

  ~FPReader () { if ( fp ) fclose (fp); }

  bool had_errors () { return error; }

private:
  FILE *fp;
  bool error;

  size_t do_fread (void *ptr, size_t size, size_t nmemb)
  {
    size_t ret = fread (ptr, size, nmemb, fp);
    if ( ret != nmemb )
        error = true;
    return ret;
  }
};

template <int N>
class differator
{
public:
  differator (const vect3d& resolution) : res(resolution) { }
  virtual ~differator () { }

  virtual void getValues (const vect3d& pos, prec_t *vals) = 0;

  void getDifference (const vect3d& pos, prec_t* dfdx, prec_t *dfdy, prec_t* dfdz,
      prec_t* d2fdxdy, prec_t* d2fdxdz, prec_t* d2fdydz,
      prec_t* d3fdxdydz)
  {
    vect3d p1, p2, p3, p4, p5, p6, p7, p8;
    prec_t h, h1, h2, h3;
    prec_t fp1[N], fp2[N], fp3[N], fp4[N], fp5[N], fp6[N], fp7[N], fp8[N];

    p1 = pos - vect3d(res.x, 0, 0);
    p2 = pos + vect3d(res.x, 0, 0);
    getValues (p1, fp1);
    getValues (p2, fp2);
    h = p2.x - p1.x;
    for ( size_t i = 0; i < N; i++ )
        dfdx[i] = (fp2[i] - fp1[i]) / h;

    p1 = pos - vect3d(0, res.y, 0);
    p2 = pos + vect3d(0, res.y, 0);
    getValues (p1, fp1);
    getValues (p2, fp2);
    h = p2.y - p1.y;
    for ( size_t i = 0; i < N; i++ )
        dfdy[i] = (fp2[i] - fp1[i]) / h;

    p1 = pos - vect3d(0, 0, res.z);
    p2 = pos + vect3d(0, 0, res.z);
    getValues (p1, fp1);
    getValues (p2, fp2);
    h = p2.z - p1.z;
    for ( size_t i = 0; i < N; i++ )
        dfdz[i] = (fp2[i] - fp1[i]) / h;

    p1 = pos + vect3d(-res.x, -res.y, 0);
    p2 = pos + vect3d(-res.x, res.y, 0);
    p3 = pos + vect3d(res.x, -res.y, 0);
    p4 = pos + vect3d(res.x, res.y, 0);
    getValues (p1, fp1);
    getValues (p2, fp2);
    getValues (p3, fp3);
    getValues (p4, fp4);
    h1 = p3.x - p1.x;
    h2 = p2.y - p1.y;
    h = h1*h2;
    for ( size_t i = 0; i < N; i++ )
        d2fdxdy[i] = ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i]))/h;

    p1 = pos + vect3d(-res.x, 0, -res.z);
    p2 = pos + vect3d(-res.x, 0, res.z);
    p3 = pos + vect3d(res.x, 0, -res.z);
    p4 = pos + vect3d(res.x, 0, res.z);
    getValues (p1, fp1);
    getValues (p2, fp2);
    getValues (p3, fp3);
    getValues (p4, fp4);
    h1 = p3.x - p1.x;
    h2 = p2.z - p1.z;
    h = h1*h2;
    for ( size_t i = 0; i < N; i++ )
        d2fdxdz[i] = ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i]))/h;

    p1 = pos + vect3d(0, -res.y, -res.z);
    p2 = pos + vect3d(0, -res.y, res.z);
    p3 = pos + vect3d(0, res.y, -res.z);
    p4 = pos + vect3d(0, res.y, res.z);
    getValues (p1, fp1);
    getValues (p2, fp2);
    getValues (p3, fp3);
    getValues (p4, fp4);
    h1 = p3.y - p1.y;
    h2 = p2.z - p1.z;
    h = h1*h2;
    for ( size_t i = 0; i < N; i++ )
        d2fdydz[i] = ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i]))/h;

    p1 = pos + vect3d(-res.x, -res.y, -res.z);
    p2 = pos + vect3d(-res.x, -res.y, res.z);
    p3 = pos + vect3d(-res.x, res.y, -res.z);
    p4 = pos + vect3d(-res.x, res.y, res.z);
    p5 = pos + vect3d(res.x, -res.y, -res.z);
    p6 = pos + vect3d(res.x, -res.y, res.z);
    p7 = pos + vect3d(res.x, res.y, -res.z);
    p8 = pos + vect3d(res.x, res.y, res.z);
    getValues (p1, fp1);
    getValues (p2, fp2);
    getValues (p3, fp3);
    getValues (p4, fp4);
    getValues (p5, fp5);
    getValues (p6, fp6);
    getValues (p7, fp7);
    getValues (p8, fp8);
    h1 = p5.x - p1.x;
    h2 = p3.y - p1.y;
    h3 = p2.z - p1.z;
    h = h1*h2*h3;
    for ( size_t i = 0; i < N; i++ )
        d3fdxdydz[i] = (((fp8[i] - fp7[i]) - (fp6[i] - fp5[i])) - ((fp4[i] - fp3[i]) - (fp2[i] - fp1[i])))/h;
  }

private:
  vect3d res;
};

#ifdef __linux__
#include <sys/time.h>

struct timer
{
  timer ()
  {
    reset();
  }

  void report (const char *name)
  {
    struct timeval tv2;
    gettimeofday (&tv2, 0);
    printf ("%s time: %u\n", name, (unsigned int)((tv2.tv_sec - tv.tv_sec) * 1000000 + tv2.tv_usec - tv.tv_usec));
    reset();
  }

  void reset ()
  {
    gettimeofday (&tv, 0);
  }

  struct timeval tv;
};
#endif

#endif // !__utils_hpp__

